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R.W. Baibieri and P.S. Schopf 
ABSTRACT 

The Kalman filter is a data-processing algorithm with 
a distinguished history in systems theory. It’s application 
to oceanographic problems is in the embryo stage and it 
is in this sense that we have demonstrated the behavior 
of the filter in the context of an internal equatorial 
Rossby wave propagation problem. 
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OCEANOGRAHIIC APPLICATIONS ON THE 
KALMAN FILTER 


INTRODUCTION 

The detennination of the state of the ocean is a major objective of descriptive physical 
oceanography, having been of concern since the earliest investigations of the seas. Even the diag- 
nosis of the mean circulation is difficult, eluding oceanographers due to an inherent variability 
in the flow and a paucity of relevant data. Limited by ship speeds and costs, the collection of 
enough in-^situ data to directly solve important circulation questions has been and will continue 
to remain as insurmountable obstacle. 

Satellite remote sensing will provide vastly increased amounts of information about the 
surface conditions of the ocean. Whether this synoptic surface data can be extrapolated to make 
quality inferences about the circulation at depth remains the fundamental question of satellite 
oceanography. Some quantities which are measured from space (e.g. SST) may be manifestations 
of deeper processes, but it is not at all obvious how one could invert the relationship to make 
reliable statements about what must be causing the surface pattern. On the other hand, there 
are some particulariy promising measurements which seem much more closely related to deeper 
flow, such as surface topography. Others, such as wind stress provide important data on the 
ocean’s forcing functions. It is felt that there is an essential link between surface topography 
and deep flow which can be exploited to make the elusive problem of ocean circulation diagnosis 
more tractable. 

The problems remaining, however, are quite formidable. The first and perhaps foremost 
difficulty is that the problem cannot be viewed as simply deterministic. It is never possible to 
attach an unambiguous, deterministic meaning to the state of the ocean from knowledge of 
previous data — initial conditjoi'.s cannot be known precisely and the equations of motion, 
however well refined, will always remain questionable, because of unknown parameters and 
un-modeled, overlooked processes. In addition, all measurements of the state of the ocean are 
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iK^, possibly biased, and suffer from sensor inaccuracies. Moreover measurements n^ade of 
the ocean are sparse in space and sparse in time so that aliasing of different length and temporal 
scales will occur. 

In this report we consider stochastic processes and observations of them. The ocean diagnosis 
problem is recast as a statistical one - given a law approximating the evolution of a stochastic 
process and given sparse, noise corrupted observations of the process, how does one construct a 
“best” estimate of the process. The word best is in quotes because the estimate changes according 
to the criteria (the performance functional) chosen to evaluate the success. 

HISTORICAL PERSPECTIVE 

The statistical tools which have been chosen to treat the oceanographic problem have their 
roots in the least-square method of Gauss in the latter part of the 18th century. He had antici- 
pated the maximum likelihood method which was presented in 1912 by RA. Fisher ( 1). Then 
Kolmogorov (2) and Wiener (3) independently derived a linear minimum mean-square estimation 
technique. This mathematical theory of smoothing and prediction was of far-reaching importance 
in communication theory. The actual computation of the optimal filter was not a simple matter 
because the weighting coefficients of the Wiener-Hopf equation had to be solved for and this in 
turn required the inversion of very large matrices. Bode and Shannon (4), simplified the theory 
and made it applicable to a wider selection of engineering problems. They realized that the Wiener- 
Hopf equation, which was essentially intractable at that time, is solvable in a fairly straightforward 
manner when a shaping filter is introduced so as to give a more explicit description of the signal 
The shaping filter yields a diHerential equation driven by white noise, the solution of which has the 
same characteristics as the signal. It was ir. 1961 that Kalman (5) presented a new approach to the 
problem of optimal filtering; his woric was based on the Bode-Shannon representation of random 
processes and the state transition techniques. One important aspect of his approach is that he 
specified the optimal estimate to be the solution of a differential equation whose coefficients are 
determined by the statistics of the process. 


2 



NEED FOR A FILTER 

It is proposed in this report that the technique which should be used to attach this estimation 
problem is the Kalman filter. First of all a filter* is needed because the dependent variables of the 
model are not always measured directly and consequently some means of inferring these variables 
from the available data must be generated. This process itself, that of inferring parameter values 
from data, is complex because the dynamic system (the ocean) is driven by forces and nonlinear 
couplings which are not well modelled and, the relationships among ^he various dependent variables 
(which describe the state of the ocean) and measured parameters are known only within some 
degrees of uncertainty. Furthermore there may be data available from many different sensors, 
each with it’s own particular dynamics and error characteristics. It would certainly be advantageous 
to have an algorithm which handled data from different sensors in an optimal way. 

There are many ways to define optimality and these are dependent upon the criteria chosen 
to evaluate performance: maximum likelihood, minimum variance, least squares, etc. For Gaussian 
noise processes the Kalman filter is optimal with respect to almost any criteria which makes 
physical sense. This will be elucidated below. 

The recursive nature of the filter means that data available at each succeeding time can be in- 
corporated with the results of the most recent filtering seep. Thus, there is no need to store all 
the oceanographic data and process it at once. 

In sununary, the Kalman filter combines all available data with a model of the dynamic system 
and prior knowledge of the physics of each sensor and produces an estimate of the desired variables 
in a way that minimizes the error of the estimate in a statistical sense. 

A discussion of the filter follows next along with a simple application for illustrative purposes. 

THE CENTRAL PROBLEM 

Irrespective of how it is formulated, the central problem of estimation theory is as follows: 
given a time varyii^ function z(t) representing a signal y(t) contaminated by noise Kt), extract 

*Frora estimation theory literature, a filter is a linear transfoimation on a stochastic process. 
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y(t) from z(t) as “cleanly” as possible when only y(t) + KO is observed. In this problem y(t), Kt) 
and their sum are considered tr be random varialles. With y(t) considered as a stochastic process, 
the probability that z(t) will attain a certain value can be computed. Furthermore given a sequence 
|z(t|)| 1 observations of the stochastic process it is in principle possible to compute the 

conditional probability distribution of y(t) given |z(t{)| j ; such a distribution represents all the 
information about y(t) contained in the sequence of measurements. This distribution yields an 
estimate y(t) of the true signal y(t); this estimate is a function of the conditional distribution and 
may be average value or the median value, peak value, etc. Up to this point the discussion has 
focussed on a reliance upon observation of the stochastic process to leam about the true signal 
y(t). But there is available another source of information about the stochastic process and that is 
dynamical models; these represent approximations to the actual dynamics of the process. The 
problem now becomes how to combine two independent sources of information about the stochastic 
process to yield an estimate of the process which is better than if only one source of information 
were used. As a preview of what is to follow, these two sources of information are combined in a 
linear, unbiased, least square sense at every discrete time point with the numerical model providing 
information not only at observation times but also in between such times. Consequently, contact 
with the stochastic process (the ocean) is never lost; at all times there is model output aiid/or 
observations available. 

Modeling of the ocean’s dynamics with analytic or numerical models is a difficult undertaking. 
Many complex physical mechanisms are either not understood, or must be eliminated from con- 
sideration for practical reasons. In addition, the initial conditions for any model calculation are not 
precisely known in a deterministic sense. 

In spite of such obstacles, analysis and estimation of ocean dynamics can still proceed, but must 
be done from a different point of view. The error sources mentioned above can be described in a 
statistical sense, and this information used in the model prediction. This is already implicitly done 
by numerical modelers. They ignore the details of sub-grid processes because they have confidence 
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that thie energy in these processes is low and that their interaction with the larger scale features is 
adequately represented through certain parameterization schemes. The statistical information 
about the error sources is passed to the system through a model noise term in the governing dif- 
ferential equations. This is usually interpreted as a white Gaussian stochastic process. (State 
augmentation techniques are readily available for handling non-white sources; the resulting aug- 
mented system model is then driven by white noise so that the procedure derived below is quite 
general and can encompass random bias noise, Markov processes, etc.) 

It will be assumed that the physical phenomena is adequately modelled by a linear stochastic 
differential equation of the form 

d = A (t) Xj dt+ B(x^,t)u(t)dt + G (x^.t) d w(t) (1) 

where A, B, C are appropriately dimensioned matrices, x^ is the stochastic process, u(t) is a deter- 
ministic forcing function and dw(t) is a Brownian motion process representing the uncertainties 
discussed above. In addition to the information about the system as provided by the numerical 
model, there are observations of the dynamical system. The observations, denoted by a vector 
z(t), are related to the stat^ of the system, x(t), by the measurement equation 

z(t) = H(t)x(t) + Kt) (2) 

where H(t) is an appropriate matrix and i^(t) is the noise embedded in the measurement. (Note 
that the transformation H allows us to predict velocity (x) and measure surface height (z), for 
instance.) The solution of equation ( 1) is not a trivial task and it is not the intent of this report 
to delve into the methods of it’s solution. For that, there are several excellent texts available; 
Gihman and Skorohad (6), Arnold (7) and, Soong (8). Suffice it to say that the solution of ( 1) 
is given by 

t t 

x(t) » *(t,tQ) x(tg) + / 4»(t,T) B(r)u(T)dr + / $(t,r) G(T)dw(r) (3) 

*o *0 
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where $ (t,t^) is the state transicion matrix. Realizing that this algorithm will be implemented on 
a digital computer, the equivalent discrete time form of (3) is given by 

x(ti^.i) » 4*(tj^j , tj) X (tj) + Bjj(tj)u(tj) + o(tj) (4) 

where it has been assumed that u(t) is constant over each interval tj<t<t^^ j and 

*i+l 

Bd(ti)= J ^tj,i,ti)B(T)dr (5) 

‘i 

*1+1 

<Xti) » <Kti+j,T)G(T)dw(T) (6) 


Furthermore a(t|) is a white Gaussian discrete time process with the same statistic as those of the 
stochastic process dw(t), namely 

E|ot(ti)|«0 (7) 

E ja(ti)a(tj)| = Q(ti)5y (8) 

In addition to the discrete time equation (4) with associated statistics (7) and (8), the measurement 
equation (2) must also have a discrete time form 


z (tj) = H(tj) x(tj) + Ktj) 

where the measurement noise vector Ktj) is assumed to have the statistics given by 

E |*<ti)| = 0 
E|Ktj) *'(tj)|= R(tj)5« 

The stage is now set for the discrete version of the Kalman filter. 


(9) 


( 10 ) 

( 11 ) 
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THE KALMAN FILTER 



The Kalman filter is a data processing algorithm which has embedded in it the dynamics of the 
system; the overall structure of the entire problem is shown in Figure 1 


DYNAMIC 

SYSTEM 

ERROR 

SOURCES 

AND 

FORCING 


SYSTEM ESTIMATE OF 



Figure 1. 


The estimation problem associated with the system equations (4) and (9) can be stated as 
follows: determine an estimate, ^(t|), of the true but unknown state x (tj) under the condition 
that this estimate should be a linear function of all observations and must satisfy the foUowing 
conditions 

(i) ^(tp must be unbiased 

(ii) il(tj) minimizes the expected value of the square of the error magnitude; 
that is, E |l(x(tj)-ll(tj)]'^ [(x(tj) - X(tj)l| » min. 

The estimation problem as just stated embodies three sub-problems: 

Filtering Problem ~ the best estimate of the state is required at the same time as 
the last measurement. 

Smoothing Problem the best estimate of the state is required at some time before 
the last measurement. 

Prediction Problem the best estimate of the state is required at some time after 
the time of the last measurement. 

The estimation problem can be viewed in terms of computations at measurement times and be- 
tween measurements. Certain variables must be propagated between measurements: the state of the 
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dynamic system and the error covariance matrix. At the time of each measurement the propagated 
state estimate must be appropriately weighted and linearly combined with the measurement to 
yield a new estimate of the state. This new estimate is then used as an initial condition for the 
dynamic system and the process begins again. 

The algorithm depends upon the estimation error defined as 

xftj) ■ ll(tp - x(tj) and x (t{) ■ 5i(tJ) - x (tj) ( 12) 

where the superscripts (-f) denote times immediately prior to and after a discrete measurement 
is incorporated at tj. Given a prior estimate of the state at time tj, what is sought is an updated 
estimate based on use of the measurement z(t|). So as to avoid a growii^ memory require- 
ments, the updated estimate is required to be a linear combination and recursive; that is the updated 
estimate must have the fonn 

8(t^»K'(t,)»(Q + K(ti)z(ti) (13) 

where K'( t|) and K(t|) are uiUcnown time varying weighting matrices. Equations (9), ( 1 2) and ( 1 3) 
are used to construct a lespiesentation for x(t|); 

x(tp *iK'(ti) + K(ti)H(tj) -I] x(ti) + K'Xtj)x(tT) + K(tj)Kt,); (14) 

equation ( 10) and the condition that the estimate i^t|) must be unbiased leads tc the representation 

K'(tj) • I - K(tj) H(tj) (15) 

With this representation for K'(t), the updated estimate (the post measurement estimate) takes the 
form 

5(tJ) • « (tJ) + K(tj) Izftj) - H(tj) ^(tpi (16) 
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The error recursion is given by 

x(tj)- [I - K(tj) H(tj)] x(t”) K(tiMti) (17) 

where K(t{) is the Kalman gain matrix. It is obtained by minimizing the sum of the diagonal 
elements of the error covariar:ce matrix 

P(t[)sEjx(t]‘)x'''(tj)j (18) 

This operation minimizes the* length of the error vector l{(t^). The resulting form for K is 


K(tj) » nQ nT(tp [H(ti)P(tpHTrtj) + R(ti)]-1 


(19) 


where superscript T denotes transpose. P(tp is the error covariance matrix prior to the incorpor* 
ation of the measurement at t| and is defined by 

P(tp « e| x(tp x”^ (tp( . (20) 


This matrix is computed during the prediction step and evolves from the last observation time by 

Vi) P(t^_i) 4»T(ti, tj.p + (Xti.i) (21) 

This best estimate of the state at t]* is obtained from 

«(tp«<Kti,ti_i)ft(t^.j) + Bd(ti_pu(ti.i) (22) 

All of the equations necessary to implement the Kalman filter have been presented, in summary* 
they are 


x(ti+j) » <^(tj+J,tJ)x(tJ) + Brt(tJ)u(tJ) + a(tJ) 
z(tt+p * H(ti+i)x(tj+i) + Kti+i) 

^(t7+p * ♦(tj+j.tpX (tJ) + Bj(tj)u(tj) 

P(t"+i) » ♦(Vi.tj)P(t|)«b'*‘(tj+j,tj) + Q(tj) 

^(t*+j) « 5(tp.i) + K(tj^.j) (z(tj+j) - H(tj+j) X(t7+j)J 
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The first two equations are state equation and the measurement equation respectively. The third 
and fourth equations represent the propagated state estimate and propagated error covariance 
respectively. The third equation may be interpreted as an estimate of the state at time tj.^j prior 
to the inclusion of an observation. -\t|.^ j ) provides a statistical measure of our confidence in the 
estimate, x(tj.^ j ), of tlie true state x( tj^ j ) prior to the use of the observation at t|^ I . It can be 
interpreted geometrically as giving the magnitude of the major and minor axes of an n-dimensional 
ellipsoid inside of which the true state x(tj.j.j) resides. The last three equations have to do with the 
time immediately after a measurement has been processed. It is seen that x(tj^. j ) is a linear com- 
bination of the propagated state and a vector z(ij^ j) - H(tj+j ) x (tj.^i ) which is known as the 
residual vector. This residual vector is the difference between the actual measurement at tj+ j , and 
the calculated measurement H(tj.j. ; ) x(tj^ j ). The matrix P(t^ j ) provides a measure of our con- 
fidence in the estimate, S(tj.j.|), of the true state x(tj.^j ) after the measurement at tj^ j has been 
incorporated. It has the same geometrical interpretation as P(t7+ j). In most instances the volume 
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of the error eUipsoid corresponding to Pltjj. j ) is smaller than that volume corresponding to P(tj+ j ) 
because the measurr ment is expected to add ^formation about the system and thereby improve 
the estimate of the state. The gain matrix K(tj^ j ) is seen to depend upon the pre-observation 
covariance matrix P(t ), and the measurement equation matrices; the gain matrix is the weight 
given to the residual vector and therefore represents the confidence to be attached to the measure- 
ments. Both P(tj^ I ) and K(t|^ I ) can be expressed in forms which reveal the effect of measurement 
noise more explicitly; that is. 

The last term in (t,> j ) is the Fisher information matrix which plays such a key role in the 
problem of observability ; it provides a measure of the contrit ition, from the observations alone, to 
the reduction in uncertainty in the estimate. The gain matrix is directly proportional to the 
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uncertainty in the estimate of the state and inversely proportional to the measurement noise. Thus 
when the measurement noise is large, confidence in estimate of the state is not significantly im- 
proved and the weight given to the measurement is small. On the other hand, when R(t^^ | ) 

+ _ 

is small, P(t ) becomes much smaller than PCtj^^) and the gain matrix is large indicating much 
more confidence in the measurement. 

A fonnal derivation of the filter equations is given in McCarty’s text (9) and is based on the ap- 
proach taken by Kalman (5), that is, from the vector space point of view and the use of orthogonal 
projections. The text by Maybeck( 10) has the derivation from the point of view of Bayes theory. 

Several aspects of what has so far been presented will now be discussed. First the statisticai be- 
haviour of the noise vectors oc(t|) and v(t^) is assumed to be white. What this means is that the noise 
is uncorrelated in time. When looked at from the power spectrum point of view, the assumption 
means that the power of the stochastic process is distributed uniformly over all frequency compon- 
ents. For the applications considered upon in this report it is known that this uniform distribution is 
not true. The assumption of white noise, although theoretically defenseless, does lead to a much 
simpler algorithm. A more meaningful approach is to let the model and/or measurement noise be 
non-white, assume a particular power spectrum (or autocorrelation) and deme a system of differ- 
ential equations corresponding to such a spectrum, which do have white noise sources. 

The procedure for handling correlated noise is presented in Maybeck ( 10) and Gelb (11). 

In equation (3) there appears the operator which takes the initial state x(t^) into the state x(t); 
this operator is known as the state transition matrix and it exists only for linear time-invariant or 
linear time-varying differential equations. In the time-invariant case 

$(t,tQ) = e^^*'*o\ 

The evaluation of the exponential as a power series may present convergence related problems if 
t-tg is not small. A closed form expression for the exponential can be derived either as a special case 
of the study of the functions of a matrix or by an algebraic method based on the Laplace transform. 
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The Laplace tran* /nn procedure is attractive when the number of equations is small. 


In the time-varying case 

and similar computational problems as in the previous case are found. 

In those car's when the non-linear equations must be used, there is available the extended 
Kalman-Bucy filter. The key aspect of the extended Alter is that the non-linear system is numeri- 
cally integrated to update the state in between measurements. 

The application of the extended filter to oceanographic problems will be the Object of our 
next paper. 

AN EXAMPLE 

To show how the Kalman Alter algorithm works, it was decided to apply it to the estimation 
of internal equatorial Rossby wave propagation. The equation of motion u*. 


9 rliv 9£v , 




A separation of variables leads to the equation 


a. 

9t 


-(2m+l)v]= -|^ 


( 22 ) 


(23) 


where m is the meridional equatorial mode. 

The problem was cast in a zonally periodic domain, with 10 grid points used to resolve the zonal 
structure. In order to keep the problem as simple as possible only two-time level marching techni- 
ques, with trapezoidal or forward timesteps, were used to integrate the model equation. With the 
second order horizontal differences, the forward scheme is unconditionally unstable. 
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The model errors are taken to be the numerical error only; no additional model noise was 

introduced as was done in Ghil, Cohn and Dalcher ( 12). The first problem was to determine the 

model error covariance matrix Q. This can be done in this case because of the luxury of having an 

analytical solution. One time step was taken and the average of the difference between the analytical 

and numerical solutions was computed as a function of the separation distance. In addition to 

having the time history of the actual errors available, denoted by A'*' and A~ in the figures, and 

+ 

P~ are also computed and these denote the traces of the covariance matrices given by P(t|) and 
P(T j) respectively. Finally each figure has a history of A™ , which is the error incurred by running 
the numerical model without incorporating any data. Figure 2 shows results using a trapezoidal 
scheme; measurement error is 2.4%, model error is 4%, initial conditions are known to within 1% 
and 5 of the 10 points are observed. Notice the history of the matrices A~, P~, these are measures 
of the actual and computed errors before data is incorporated. There errors are stable and hover 
around the 5% error level. Now notice the error matrices A”*', P"*", which are measures of the actual 
and computed errors after measurements have been incorporated. These erron are also stable and 
hover around the 2% level. It is clear how the use of measurements serves to reduce the pre- 
measurement covariance metrics, A~, P~, and how the application of the Alter serves to significantly 
reduce the error made by leaving the numerical model alone and not applying the Kalman filter at 
aU. 

Figure 3 shows results when only the initial condition uncertainty was changed. Here, the 
initial conditions are not known very well with an error of about 3 1.6%, compared to 1% in Figure 
2. The results the same as in Figure 2 ; here it is clearly shown how the filter behaves when 
initial conditions are poorly known. This behavior is attributed to the fact that the measurement 
noise is small. The initial error has been erased after two observations. 

Figure 4 shows results when only one grid point has been observed. The input conditions are 
otherwise the same as in Figure 2. Several points are clear; first the excursions of A~ and A'*’ about 
P~ and P'*' are larger than in the previous two figures. Secondly the incorporation of information 
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from one-grid point hardly reduces the uncertainty P'*' of the estimate of the state. The filter 
results P“, P"*" asymptotically approach a level of about 7%. 

Figure S shows what happens when only one point is observed and the initial condition error 
is large. As in Figure 3, the filter takes slightly longer to erase the large uncertainty in initial con- 
ditions; in particular five measurements at the one grid point are sufficient to bring the filter per- 
formance down to the same error level as in Figure 4. 

Figure 6 shows results when the forward timestep scheme is used. The numerical model be- 
comes unstable, and consequently the estimate of the model noise error covariance matrix, Q, is 
slightly larger than in the first four figures; in particular the model noise is estimated after one 
time step to be 7.2%. All other input parameters are the same as in Figure 2. The results show 
that both A“ and P~ track very well with an rms error level of about 8%, as opposed to about 
5%, in Figure 2. This level of error is only slightly larger than the assumed model error. When data 
is incorporated the rms level of uncertainty of the state estimate approaches 3%. However there is 
a hint that after many time steps the post-measurement actual error, as given by A'*’, begins to 
diverge from the post-measurement error as computed by the filter and given by P"*". This is an 
encouraging finding for it says that even though the numerical integration is unstable and will 
diverge when no data is used (see the plot of A*"), the incorporation of data and the Kalman 
filter prevents the unstable scheme from getting away. The use of data yields a new estimate of 
the state of the dynamic system and this new state is used as an initial condition for the next step. 
It is the constant updating of initial conditions that brings the results of the unstable scheme back 
so that at the beginning of the next step the numerical model starts closer to the analytical model. 

Figure 7 shows what happens when the forward timestep scheme is used with a large initial 
condition error and only one point observed. During the first several time steps the filter is able to 
significantly reduce the initial conditions error. However the observation of only one point is not 
sufficient to prevent the un?^..ble scheme from influencing the overall filter performance and con- 
sequently, the error covariance matrices become very large. 
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CONCLUSIONS 

The extraction of information from measurements of a dynamic system governed by a 
stochastic evolution law is a characterization of the general ocean circulation problem. The Kalman 
filter, a specific tool from the framework of recursive estimation theory, provides an optimal way 
of utilizing measurements and model output to obtain greater confidence in the estimate of the 
systems evolution. 

The filter has been applied to the internal equatorial Rossby wave propagation problem with 
rather encouraging results. Analysis of numerical experiments has demonstrated the filter’s per- 
fonnance in the presence of initial condition errors, measurements errors and modelling errors. 

For the class of problems treated here, initial condition errors are rapidly erased, within 5 time 
steps. Measurement errors tend to add a constant bias to the covariance matrices. 

The most encouraging result is the performance of the filter in the presence of a numerically 
unstable integration scheme. Provided there are enough observations of the dynamic system, the 
filter tends to weight these observations more heavily and the unstable scheme never has the 
opportunity to overwhelm the filter and force divergence. 

Also demonstrated is the performance of the filter using different observation schemes. When 
only one grid point is observed, the actual errors tend to show slightly larger deviations than the 
case when five grid points are observed. The covariance matrices P“ and P"*" show the same shape 
time history as the case when five points are observed; however the asymptotic value is slightly 
larger. In summary, the Kalman filter and its extensions show much promise as a tool for analyzing 
oceanographic problems from a statistical point of view. 
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